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13.0  FINITE  ELEMENT  FORMULATION  OF  THE  LAYER-WISE  SHEAR  DEFORM- 
ABLE  COMPOSITE  SHELL  THEORY 

13.1  Introduction 

Present  Composite  Shell  Technology  and  its  Limitations 

13.1.1  Qassical  lamination  theories 

Plate  and  Shell  structures  made  of  laminated  composite  materials  are  often  noodeled  as  an 
equivalent  single  layer  using  classical  laminate  theory  (C.L.T.)  in  which  the  thickness  stress 
components  are  ignored.  The  classical  laminate  theory  is  a  direct  extension  of  classical  plate 
and  shell  theory  in  which  the  well  known  Kirchhoff-Love  kinematic  hypothesis  is  enforced. 
This  theory  is  adequate  when  the  thickness  (relative  to  side  or  radius)  is  small  so  that  the  varia¬ 
tion  of  the  field  variables  through  the  thickness  direction  is  minimal.  However,  laminated  plates 
and  shells  made  of  advanced  filamentry  composite  materials  are  succeptible  to  thickness  effects 
because  their  effective  transverse  moduli  are  significantly  smaller  than  the  effective  elastic 
noodulus  along  the  fiber  direction.  Furthermore,  the  classical  theory  of  plates  which  assumes  that 
the  normals  to  the  midplane  before  deformation  remain  straight  and  normal  to  the  plane  after 
deformation,  underpredicts  deflections  and  overpredicts  natural  firequencies  and  buckling  loads. 
These  discrepancies  are  due  to  the  neglect  of  transverse  shear  strains.  The  errors  in  deflection, 
stresses,  natural  frequencies,  and  buckling  loads  are  even  higher  for  plates  made  of  advanced 
composite.  TTie  range  of  applicability  of  the  C.L.T.  solution  has  been  well  established  for  lam¬ 
inated  flat  plates  by  Pagano  [see  Pagano  1989].  These  analyses  have  indicated  that  a  theory 
which  accounts  for  the  transverse  shear  deformation  effects  would  be  adequate  to  predict  the 
gross  behavior  of  the  laminate. 

13.1.2  Shear  deformation  theories 

In  order  to  overcome  these  deficiencies  in  C.L.T.,  refined  laminate  theories  have  been  pro¬ 
posed.  These  are  single  layw  theories  in  which  the  transverse  shear  stresses  are  taken  into 
account  They  provide  improved  global  response  estimates  for  deflections,  vibration  fiequencies 
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and  buckling  loads  of  naoderately  thick  composites  when  compared  to  the  classical  laminate 
theory.  A  Mindlin  type  first-order  transverse  shear  deformation  theory  (S.D.T.)  was  first 
developed  by  Whimey  and  Pagano  [1970]  for  multilayered  anisotropic  plates,  and  by  Dong  and 
Tso  [1972]  for  multilayered  anisotropic  shells.  Both  of  these  approaches  (CX.T.  and  S.D.T.) 
considered  all  layers  as  one  equivalent  single  anisotropic  layer,  thus  these  approaches  are  inade¬ 
quate  to  model  the  warpage  of  cross-sections,  that  is,  the  distortion  of  the  deformed  normal  due 
to  transvase  shear  stresses.  Furthermore,  the  assumption  of  nondeformable  normal  results  in 
incompatible  shearing  stresses  between  every  two  adjacent  layers.  Also  the  latter  approach 
requires  the  introduction  of  an  arbitrary  shear  correction  factor  which  is  dependent  on  the  lami¬ 
nation  parameters  for  obtaining  accurate  results. 

13.1.3  3-D  Anisotropic  elasticity 

In  another  class  of  composite  problems,  i.e.,  for  composite  structures  with  a  thick  cross- 
section,  two-dimensional  plate  analyses  are  inadequate  because  through-thickness  stresses  (inter¬ 
laminar  and  normal  stresses)  are  comparable  in  magnitude  to  the  other  stress  components.  Thus 
a  three-dimensional  finite  element  analysis  is  necessary  in  order  to  calculate  the  through¬ 
thickness  stresses  accurately.  Since  material  properties  vary  from  layer  to  layer  due  to  the 
change  of  the  ply  orientation,  finite  element  modeling  for  very  thick  composites  throughout  the 
thickness  becomes  extremely  difficult  and  expensive.  Traditional,  three-dimensional  finite  ele¬ 
ment  methods,  based  on  one  layer  per  element,  are  not  computationally  efficient  for  analysis. 

New  Directions  in  Composite  Shell  Technology 

The  exact  analyses  performed  by  Pagano  [1989]  on  the  composite  flat  plates  have  indicated 
that  the  distortion  of  the  deformed  normal  is  dependent  not  only  on  the  laminate  thickness,  but 
also  on  the  orientation  and  the  degree  of  orthotropy  of  the  individual  layers.  Therefore  the 
hypothesis  of  nondeformable  normals,  while  acceptable  for  isotropic  plates  and  shells  is  often 
quite  unaccqrtable  for  multilayered  anisotropic  plates  and  shells  with  very  large  ratio  of  Young’s 
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modulus  to  shear  modulus,  even  if  they  are  relatively  thin.  Thus  a  transverse  shear  deformation 
theory  which  also  accounts  for  the  distortion  of  the  deformed  normal  is  required  for  accurate 
prediction  of  the  behavior  (deflections,  thickness  distribution  of  the  in-plane  displacements, 
natural  frequencies,  etc.)  of  multilayered  anisotropic  plates  and  sheUs. 

In  view  of  these  issues  a  variadonally  sound  theory  that 

•  Accounts  for  the  3-D  effects 

•  Allows  thickness  variation,  and 

•  Permits  the  warping  of  the  deformed  normal 

is  required  for  refined  and  sophisticated  analysis  of  thick  and  thin  composites. 

The  approach  proposed  in  this  work  utilizes  a  displacement  field  which  fulfills  a  priori  the 
static  and  geometric  continuity  conditions  between  contiguous  layers.  It  is  worth  mentioning  that 
the  number  of  partial  differential  equations  in  the  resulting  system  is  independent  of  the  number 
of  plies.  In  addition,  the  order  of  the  system  is  the  same  as  in  the  first-order  shear  deformation 
theory.  The  chief  advantage  of  the  assumed  displacement  field  rests  on  its  capability  to  model 
the  distortion  of  the  deformed  normal  and  to  satisfy  the  contact  conditions  ab  initio,  without 
increasing  the  number  and  order  of  the  partial  differential  equations  with  respect  to  the  first- 
OTdCT  transverse  shear  deformation  theory.  Furthermore,  it  is  feasible  to  employ  this  formulation 
for  constructing  plate  and  shell  finite  elements  via  the  finite  element  displacement  method. 
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13.2  A  C^’-continuous  Theory  for  Composite  Laminates 

The  Continuum  Theory  for  Composite  Laminates  as  presented  in  the  previous  chapters 
requires  C^-continuity  of  shape  functions,  just  as  does  the  classical  Poisson-Kirchhoff  Plate 
Theory  and  the  classical  BemouUi-Euler  Beam  Theory. 

From  a  mathematical  viewpoint,  if  the  stiffness  integrands  involve  derivatives  of  order  m, 
then  the  requirements  for  the  convergence  of  the  Finite  Element  Solution  to  the  exact  solution 
with  the  refinement  of  the  mesh  are 

i)  Shape  functions  should  be  smooth  to  order  C”  on  each  element  interior,  12*; 

ii)  Shape  functions  should  be  continuous  across  each  element  boundary  P. 

iii)  Shape  functions  should  be  complete. 

Finite  elements  that  satisfy  these  properties  are  called  CONFORMING,  or  COMPATIBLE 
elements. 

Continuous  (i.e.,  finite  element  interpolations  are  easily  constructed.  The  same  cannot 
be  said  for  multidimensional  C^-interpolations  [Hughes,  87].  Furthermore,  there  has  been  an 
increasing  trend  in  the  literature  to  go  towards  elements  based  upon  theories  which  accommo¬ 
date  transverse  shear  strains  and  require  only  C°-continuity  [Hughes,  87].  Although  this 
approach  is  not  without  its  own  inherent  difficulties,  it  opens  the  way  to  a  greater  variety  of 
interpolatory  schemes. 

13.3  New  Ideas  Proposed  in  the  Present  Theory 

In  the  following  section  we  present  a  Finite  Element  Shear  Deformable  Theory  for  thick  as 
well  as  thin  composites.  The  proposed  theory  has  a  number  of  noteworthy  attributes  which 
directly  address  the  technical  drawbacks  present  in  most  of  the  theories  that  have  been  proposed 
for  composite  analysis  to  date. 
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•  The  displacement  field  proposed  in  this  work  is  continuous  in  3-D  where  as  the  rota¬ 
tion  field  is  layer  '.ise  continuous  (in  2-D)  and  can  be  discontinuous  across  the  finite 
element  laj  .^s  dux>ugh  the  thickness  direction. 

•  The  displacement  field  fulfills  a  priori  the  static  and  geometric  continuity  conditions 
between  ctmtiguous  layers. 

•  The  novel  idea  in  the  assumed  displacement  field  lies  in  its  capability  to  model  the 
distortion  of  the  deformed  normal,  without  increasing  the  number  and  order  of  the  par¬ 
tial  differential  equations  with  respect  to  the  first-order  transverse  shear  defonnation 
theory. 

•  Another  new  idea  in  the  theory  is  its  3-D  feature,  thereby  modeling  the  interlaminar 
conditions  and  predicting  the  3-D  edge  effects  more  accurately. 

•  A  salient  feature  of  the  proposed  theory  is  that,  at  most,  only  first  derivatives  of  dis¬ 
placement  and  rotation  fields  appear  in  the  variational  equations.  The  practical 
consequence  of  this  fact  is  that  only  continuity  of  finite  element  functions  is 
required  which  is  readily  satisfied  by  the  family  of  Lagrange  elements. 

•  The  number  of  partial  differential  equations  in  the  resulting  system  is  independent  of 
the  number  of  plies  and  their  orientations  in  the  composite. 

•  Another  advantage  of  the  proposed  composite  shell  theory  lies  in  the  greater  flexibility 
in  the  specification  of  the  boundary  conditions. 

•  The  theory  covers  a  wide  range  in  the  sense  that  in  one  limit  case  when  there  is  only 
one  layer  of  proposed  elements  through  the  thickness,  one  recovers  the  features  of  the 
standard  Shear  Deformation  Theories  (S.D.T.).  However  the  added  advantage  in  the 
present  case  lies  in  the  3-D  feature  of  the  theory  which  controls  the  variation  in  the 
thickness  via  the  Poisson  terms  rather  than  ad  hoc  mathematical  tricks  as  done  in  the 
literature. 

•  In  another  limit  case,  one  can  model  the  composite  with  one  element  per  ply  through 
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the  composite  thickness,  a  procedure  that  is  typically  done  while  using  the  standard 
3-D  anisotropic  elasticity  elements.  The  added  advantage  of  the  proposed  theory  in 
this  limit  case  is  that  because  of  the  shear  deformation  capability  of  the  proposed  ele¬ 
ments,  they  model  the  warping  of  the  deformed  normal  more  accurately,  thereby 
improving  the  bending  behavior. 

•  From  a  practical  design  point  of  view  it  provides  the  engineer  the  freedom  to  deter¬ 
mine  the  precision  in  analysis.  If  a  general  response  of  the  composite  structure  is 
required,  the  composite  can  be  modeled  with  one  element  through  the  thickness.  On 
the  other  hand,  the  designer  can  model  the  thickness  with  as  many  layers  of  the  pro¬ 
posed  element  as  deemed  necessary  to  achieve  the  required  accuracy. 

•  Furthermore,  it  is  feasible  to  employ  this  formulation  for  constructing  plate  and  shell 
finite  elements  via  the  finite  element  displacement  method. 

13.4  Main  Assumptions  of  the  Piate/Shell  Theories 
1.  The  domain  O  is  of  the  following  special  form; 

n  =  {(x,y,z)  €  I  z  e  [-  Y  ,  y]  ,  (x,y)  €  A  c  F?}  (13.4.1) 

where  t  is  the  plate  or  the  shell  thickness  and  A  is  the  area  of  the  reference  surface. 


2.  033  =  0 

=>  plane  stress  hypothesis  (see  Hughes,  87,  p.  31 1). 

3.  U„(x,y,z)  =  -z0a(x,y) 

This  assumption  implies  that  plane  sections  remain  plane.  6^  is  interpreted  as  the  rotation 
of  a  fiber  initially  normal  to  the  plate  reference  surface. 
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4.  U3(x,y,z)  =  W(x,y) 

This  means  that  the  transverse  displacement  w  does  not  vary  through  the  thickness. 

13J5  Main  Assumptions  of  the  Layer-wise  Sh^r  Deformable  Shell  Theory 

1.  The  domain  £2  is  of  the  following  special  form: 

(0 

,  T  =  £ z('), (x.y)<'> €  cl^ -  (13.5.1) 

t 

i.e.,  there  are  /  layers  of  finite  elements  in  the  thickness  direction,  and  is  the  area  of  the 
reference  surface  for  that  layer  and  T  is  the  total  thickness  T  of  the  composite  shell. 

2.  The  displacement  field  is  assumed  to  take  the  following  form 

U£^(x,y,z)  *  u4'^(x,y,z)  -  z%Q^(\,y)  (13.5.2) 

here  u4^(x,y,z)  are  the  displacements  and  04^(x,y)  are  the  fiber  rotations  for  the  /*  layer 
reference  surface  and  ^  6  [0,1]  is  a  parameter  that  establishes  the  position  of  a  point  from 
the  reference  surface  in  the  thickness  direction. 

3.  The  displacement  field  in  the  thickness  direction  is  assumed  to  be  a  function  of  z 

U^'^  =  U3(x,y,z)  (13.5.3) 

By  this  assumption  the  transverse  displacement,  U3  does  vary  through  the  thickness  thereby 
producing  through  the  thickness  strains  which  result  in  thickness  variation  in  the  shell. 

4.  As  a  consequence  of  the  above  relaxation 

O33  ^  0 


n=-  (x,y. 
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i.e.,  we  do  not  invoke  the  plane  stress  hypothesis. 

S.  We  know  from  the  elasticity  theory  that  the  displacements  and  stresses  at  the  interface 
between  /*  and  (/+!)*  bounded  layers  must  satisfy  the  following  contact  conditions 

=  (13.5.4a) 

=  (13.5.4b) 

In  addition,  continuity  of  stress  tractions  requires  that  the  foUowing  stress  conditions  at  the 
interface  of  layers  be  satisfies’ 


=  (I3.5.5a) 

=  (13.5.5b) 

A  consequence  of  the  second  assumption  above,  (13.5.2),  is  that  each  finite  element  layer  is 
associated  with  non-normal  cross-sectional  rotations  which  are  assumed  to  be  the  same  for  all 
representative  elements  or  plies  in  the  finite  element  layer  in  accordance  with  the  Mindlin 
kinematic  assumption.  Another  consequence  of  the  second  assumption  is  that  it  results  in 
independent  shear  deformation  of  the  director  in  each  layer  and  allows  the  warping  of  the  com¬ 
posite  cross-section.  It  also  results  in  discontinuous  strain  fields  across  the  different  material 
sets,  thereby  creating  the  provision  of  stress  continuity  across  the  material  interfaces.  Conse¬ 
quently,  tire  fifth  assumption  is  inherently  satisfied. 

The  proposed  theory  goes  one  step  further  than  the  present  shear  deformable  theories  for 
composite  laminates  in  that  the  non-normal  cross-sectional  rotations  may  or  may  not  be  the  sariK 
from  one  finite  element  layer  to  the  other.  This  allows  warping  or  deformation  of  the  normal, 
thus  producing  a  higher  order  displacement  field  through  the  thickness  direction.  Consider  the 
interface  between  and  (/+1)^  finite  element  layers  and  impose  the  continuity  of  the  tangential 


components  of  the  displacement  field  via 


u(/)=ur' 

therefore 

\x^{x,y^)  +  =  u^\x,y,z)  +  K^^^\x,y)  (1 3.5.6) 

for  the  interface  between  /  and  (/+!)*''  layer 

C  =  +1  for  /*  layer 
4  =  0  for  (/+!)*  layer 

Substituting  in  Eq.  (13.5.6) 


}x^^\x,y,z)  -  u^^(x.y,z)  =  z^%^(x,y) 

ei'^Cx.y)  -  ~  [uf ‘>(x,y.z)  -  u«>(x,y,z)]  (13.5.7) 

It  implies  that  in  a  discrete  sense 


ei%,y)  =  uS(x,y,z)  (13.5.8a) 

0iVy)  =  uS(x.y^)  (13.5.8b) 

So  when  we  substitute  these  in  Yq3  expressions  we  obtain  the  modified  shear  strain  expression. 

Furthermore,  in  examining  and  in  assumptions  2  and  3,  it  is  easy  to  realize  that  the 
displacement  field  assumed  is  a  continuous  function  of  z  coordinate  for  all  values  of  u^'^(x,y,z) 
and  0^^(x,y).  Consequently,  the  requirement  of  continuity  of  the  displacement  field  is  automati¬ 
cally  satisfied. 


0 


13.6  Geometric  Representation 
13.6.1  Proposed  Composite  Shell  Element 

Figure  13.1  shows  a  typical  configuration  of  a  doubly  curved  composite  shell.  It  can  be 
made  of  numerous  plies  with  variable  material  properties,  reinforcement  fiber  orientations  and 
ply  thicknesses.  For  all  practical  purposes,  these  plies  are  stacked  on  top  of  each  other  in  a  cer¬ 
tain  predefined  sequence.  We  call  this  sequence  which  repeats  itself  in  the  thickness  direction  as 
a  "representative  element"  In  our  mathematical  nsodeling  of  the  mechanics  of  shell,  no  such 
assumption  has  been  made  which  limits  the  number  of  individual  plies,  ply  thicknesses,  their 
orientation  or  their  stacking  sequence  to  join  a  representative  element  At  the  same  time,  there  is 
no  limitation  on  the  number  of  such  representative  elements  in  the  thickness  of  the  composite 
shell.  Consequently  the  finite  element  model  developed  can  be  applied,  with  the  same  ease,  both 
at  the  individual  laminate  level  (i.e.,  the  microstructurc  level)  or  the  representative  element  level 
(i.e.,  the  macrostructurc  level)  or  even  to  a  congregate  of  representative  elements. 
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As  shown  in  Fig.  13.1,  the  composite  is  discretized  via  four  finite  element  layers  through 
the  thickness.  As  noted  previously,  each  of  these  layers  could  be  representing  either  the  micro  or 


Fig.  13.1b 

Stress  representation  on  a  typical  finite  element  layer 


Let  us  concentrate  our  attention  to  the  case  in  which  the  composite  is  composed  of  four 
laminates,  each  being  modeled  via  a  finite  element  layer.  A  typical  such  layer  is  shown  via  the 
shaded  region  in  Fig.  13.1a.  The  stresses  and  stress  couples  acting  on  this  layer  are  shown  in 
Fig.  13.1b.  Figure  13.2  shows  a  schematic  diagram  of  the  geometry  of  the  layer  with  the  refer¬ 
ence  surface  associated  with  the  bottom  surface  of  the  layer.  In  our  presentation,  a  typical  com¬ 
posite^  shell  element  has  its  reference  surface  associated  with  its  lower  face  as  shown  in  Fig. 
13.2. 
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13.6.2  Kinematics  of  the  PefOTmarion  Through  the  Thickness 

Figure  13.4  shows  a  thick  composite  laminare  with  "m"  number  of  representative  elements, 
each  ccHitaining  "n"  different  layers.  In  the  "Layerwise  Shear  Deformable  Finite  Element 
Theory  for  Coiiq;)Osite  Shells"  that  we  are  presenting  in  this  work,  the  tmal  thickness  T  is  divided 
into  L  finite  elonents  through  the  thickness.  Each  layer  of  finite  elements  through  the  thickness 
is  associated  with  a  reference  surface  which  is  coincident  with  the  lower  surface  of  that  layer. 
The  total  thickness  is  given  by 

ti  +  t2+  •••  +t(;,=:T  (13.6.1) 

and  reference  surface  of  1st  layer  lies  at 

t  =  0 

Reference  surface  for  the  second  layer  of  finite  elements  lies  at 


t  =  ti 

Similaily  for  the  3rd  layer 

t  =  (t]+t2) 

It  is  very  important  to  note  that  the  normal  fiber  rotation  (i.e.  Gq)  and  die  slope  (i.e.,  U3  „) 
are  not  necessarily  the  same  and  thus  transverse  shear  strains  are  accommodated.  This  is  to  be 
contrasted  with  the  classical  lamination  theories  (C  J^.T.)  in  which  Gq  =  consequently  the 

transverse  shear  strains  are  zero. 

Within  each  layer  (of  finite  elements),  the  ncnmal  to  the  reference  surface  for  that  layer, 
rotates  by  Gq,  thereby  generating  shear  strain  703.  Consequently,  in  the  deformed  configuration, 
node  2  (through  the  thickness,  see  fig.  13.4)  moves  to  the  location  2'. 


Now,  the  nomtal  to  the  undeformed  reference  surface  in  the  second  layer  rotates  by  an 
angle  6q,  generating  shear  strain  Ya3-  (Here  1  =  2,  i.e.,  finite  element  layer  #2).  Because  of  the 
continuity  of  the  displacement  field  in  the  thickness  direction  we  obtain  the  new  locations  of 
finite  element  nodal  points  as  T,  2\  3'  and  4',  shown  in  the  deformed  configuration.  It  should  be 
noted  that  this  new  location  of  points  produces  a  higher  order  variation  of  strains  through  the 
thickness  direction.  This  is  a  very  important  feature  that  precludes  the  need  for  introducing  ad 
hoc  polynomial  expressions  to  model  the  higher  order  variation  of  displacement  field  through  the 
thickness. 
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13.7  Finite  Element  Description  of  the  Thick  Composite  Shell 
13.7.1  Doubly  curved  c  ^.troosite  shells  in  3-D 

The  geometry  of  a  typical  quadrilateral  shell  element  is  defined  by  the  following  relations 


"“’ft.’i.O = •"’(in) + 

(13.7.1) 

(13.7.2) 

x%n.C)=  S 

(13.7.3) 

X<'>(0  =  Z‘'>(0Xf  ("osmi) 

(13.7.4) 

z<'>(C)=N*(Ozi®'+N-(Ozi'>- 

(13.7.5) 

N*®  =  |(l+0 

(13.7.6a) 

N-ffi  =  1  (1-C) 

(13.7.6b) 

where 

the  position  vector  of  a  generic  point  of  the  shell  for  layer  /. 
x^^:  the  position  vector  of  a  point  in  the  reference  surface,  for  layer  /. 

the  portion  vector  of  a  generic  point  relative  to  x^'^  which  defines  the  director  through 
the  point  for  layer  /  (in  conq)utational  shell  literature,  is  referred  to  as  fiber  direc¬ 
tion). 

x^^:  the  position  vector  of  nodal  point  "a"  in  layer  /. 

N^:  denotes  a  two-dimensional  shape  fimetion  associated  with  node  "a". 

Den^  the  number  of  element  nodes  in  the  reference  surface  of  layer  /. 
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X^:  a  unit  vector  emanating  from  node  "a"  in  the  director  direction. 

a  "thickness  function"  associated  with  node  "a",  which  is  defined  by  the  location  of  the 
reference  surface. 

The  above  relations  represent  a  smooth  mapping  of  the  biunit  cube  into  the  physical  shell 
domain.  For  ”C"  fixed,  the  surface  defined  by  (13.7.1)  is  caUed  a  lamina  and  for  fixed,  the 
line  described  by  (13.7.1)  is  the  director.  The  directors  are,  in  general,  not  perpendicular  to  the 
laminae.  Sometimes  the  director  is  referred  to  as  the  "pseudonormal." 

For  a  particular  choice  of  two-dimensional  shape  functions,  eqs.  (13.7.1)-(13.7.6)  are  pre¬ 
cisely  defined  upon  specification  of  X^,  and  (a  =  l,2,...ne^.  It  is  convenient  in 
practice  to  take  as  input  the  coordinates  of  the  top  and  bottom  surfaces  of  the  shell  along  each 
nodal  director  (x^^  and  x^^,  respectively)  and  a  parameter  ^  6  [-1,+1],  which  defines  the  loca¬ 
tion  of  the  reference  surface.  For  example  if  4  -1,0,+1  (respectively),  then  the  reference  sur¬ 

face  is  taken  to  be  the  bottom,  middle,  top  (respectively)  of  the  shell.  From  these  data  we  may 
calculate 


=  j  y  (1+C)*i^ 

(13.7.7) 

Y(n_  *•  *• 

•  llxi^^-xi^ll 

(13.7.8) 

=  ^  (K)  llxi^-xi'>-|l 

(13.7.9a) 

^  (1+0  »  xi'^^  -  x<'>- II 

(13.7.9b) 

whae  II  •  II  denotes  the  Euclidean  norm  (i.e.,  II x II  =  (xj+x^+x^y^). 
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Lamina  Coordinate  Systems 


At  each  integration  point  in  the  element  a  Cartesian  reference  frame  is  erected  so  that  two 
axes  are  tangent  to  the  lamina  through  the  point.  The  frame  is  defined  by  its  orthononnal  basis 
vectors  in  which  ej'  is  perpendicular  to  the  lamina.  The  basis  vectors  are  calculated  as 

follows 


Rg.  13.5  Typical  lamina  coordinate  system  =  constant  surface) 
Construct  unit  tangent  vectors  to  the  and  il-  coordinate  directions: 


Consequently 


(13.7.10a) 


(13.7.10b) 


(13.7.11) 


The  vectors  tangent  to  the  lamina  are  selected  so  that  the  angle  between  ef  and  is  the  same 
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as  the  angle  between  and  ej*  and  so  that  the  eiKe^-basis  is  as  "close"  as  possible  to  the 
basis.  Thus 


=  (13.7.12a) 

=  +  (13.7.12b) 


where 


(/)__£ _ 

Ily(ef  +  e5'^)ll 


ej-xe<^ 

llej'xe^^ll 


(13.7.13a) 


(13.7.13b) 


Furthermore,  we  also  define  the  orthogonal  matrix  k>  transform  quantities  from  the  global  coor¬ 
dinate  system  to  the  lamina  system 


q  =  [qij)  =  [ei^el'ej'f 


(13.7.14) 
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For  the  present  case 


Fig.  13.6  11=  constant  surface 


Ci'ssAj  ,  e^-A2  ,  63*  =  D 
When  specialized  to  flat  geometries  these  become 


1 

=  =Ei  =i 

[  1  1 
0  > 
L  oj 

=  A2  =  E2='j 

! 

f  01 
1  ' 

1  oJ 

f 

■  01 

:L  =  d  =  E,=-^ 


0 

1 


(13.7.15) 


(13.7.16a) 


(13.7.16b) 


(13.7.16c) 


where  Ej,  E2,  E3  are  unit  vectors  in  the  Cartesian  frame. 
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Fiber  Coordinate  System 

A  unique  coordinate  system  is  erected  at  the  nodes  on  the  reference  surface  for  each  layer 
and  is  used  as  a  reference  frame  for  rotations.  One  of  the  directions  of  this  frame  is  required  to 
be  coincident  with  the  director  direction.  This  one  condition  is  not  sufficient  to  define  the  frame. 
The  following  algorithm  [Hughes  1987]  can  be  employed  to  define  this  frame. 

A 

Let  X  denote  the  unit  basis  vector  in  the  director  direction  and  61,62,03  denote  the  global 
Canesian  basis,  i.e.. 


1 

\  01 

«!=■ 

0 

L  oj 

II 

A 

_A _ 

L  Jj 

'  ♦63=- 

0 

L  1 J 

The  global  Canesian  components  of  X  are  denoted  by  Xj,  i  =  1,2,3. 

Al£orithm 

1.  Utajs  1^1,  i  =  1,2,3. 

2  j  =  l 

3.  Ifai  >33,  thena3  =  ai,  andj=2. 

4.  If  a2  >  a3,  j  =  3, 

5.  e|=X 

6.  6|=(Xxej)/IIXxejll 

7.  ef  =  e|xX 

The  orthonormal  fiber  basis  obtained  (i.e.,  ei,e2,63)  satisfies  the  condition  that  if  X  is  close 
to  63,  then  61,62,63  will  be  close  to  61,62,63,  respectively. 
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13.7.2  Kinematics  in  the  context  of  Finite  Element  Method 

The  displacement  of  the  shell  is  assumed  to  take  the  following  form; 


= ^^1]) + 

(13.7.17) 

u^'^(4.Tl)=XN.a,Ti)u« 

(13.7.18) 

S=1 

(13.7.19) 

U«(C)  =  ZP(C)U«  (no  sum) 

(13.7.20) 

where 

is  the  displacement  of  a  generic  point  in  the  shell  layer  /. 

is  the  displacement  of  a  point  on  the  reference  surface  of  the  shell  layer  /. 

UW 

is  the  "director  displacement"  for  the  shell  layer  /. 

U(/)  =  U(/)  +  u(0 

where 

The  vector  is  constructed  such  that  the  director  may  rotate,  viz. 

u<'>=e^eiP-eMf 

The  quantities  6^P  and  0^  represent  the  rotations  of  the  fiber  about  the  base  vectors  and 
e^,  for  shell  layer  /,  respectively. 
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13.7.3  Appearance  of  3D  effects  in  the  theory 

Typical  shell  theories  have  shell  description  and  parameterization  done  via  a  2D  reference 
surface  and  a  director  field  which  is  orthogonal  to  this  surface.  Plane  stress  hypothesis  is  nor¬ 
mally  invoked.  Based  on  the  assumption  that  the  thickness  of  the  shell  is  small,  some  of  these 
theories  do  not  permit  the  thickness  variation  thereby  making  the  thickness  strains  identically 
equal  to  zero.  In  other  cases.  [Simo  and  Fox  1988].  the  variation  of  the  shell  thickness  is 
accounted  for  by  associating  a  parameter  with  the  direaor  field  and  bounding  it  so  that  the  thick¬ 
ness  does  not  take  unrealistic  values.  First  allowing  and  then  later  controlling  the  thickness  vari¬ 
ation  of  shells  is  an  important  issue  and  it  arises  because  of  the  introduction  of  engineering 
approximates  to  convert  a  3D  theory,  representing  3D  phenomenon,  to  a  2D  theory  and  still 
expecting  it  to  somehow  manifest  the  effects  in  die  third  dimension. 

In  our  description  of  the  shell,  there  is  a  third  dimension.  i.e..  through  the  thickness  dimen¬ 
sion  which  permits  the  development  of  strain  field,  and  thereby  modeling  thick  laminates  in  a 
more  realistic  way.  Furthermore,  we  do  not  have  to  introduce  auxiliary  functions  in  the  director 
field  to  control  the  shell  thickness.  In  our  theory,  because  of  the  Poisson  effects  in  the  third 
dimension,  strains  develop  through  the  thickness  and  account  for  the  thickness  variation. 
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13.8  Derivation  of  C  matrix  for  a  representative  element 


Fig.  13.7  Schematic  diagram  of  a  "representative  element"  and  its  constituent  laminates 
Suppose,  the  representative  element  is  made  of  two  laminates  of  the  same  material  but  with 
reinforcement  fiber  orientation  at  +6  and  -0  degrees  with  respect  to  the  extension  direction. 

Let  C((x)  (a  =  1... .number  of  plies  in  representative  element)  represent  the  C  matrix  (consti¬ 
tutive  matrix)  for  the  laminate  with  regard  to  its  mutually  perpendicular  planes  of  elastic  sym¬ 
metry.  In  general  it  can  be  written  as 
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C(o)  = 


Cii 

C12 

Cl3 

0 

0 

^16 

C12 

C22 

C23 

0 

0 

C26 

Ci3 

C23 

C33 

0 

0 

C36 

0 

0 

0 

C44 

C45 

0 

0 

0 

0 

C54 

C55 

0 

C16 

C26 

C36 

0 

0 

^66 

(13.8.1) 


where  a  represents  the  number  of  laminates  in  a  representative  element. 


Let  represent  the  direction  along  the  loading  for  the  composite  element  and  let  X3 
represent  its  thickness  direction.  This  axis  is  assumed  to  be  perpendicular  to  the  plane  of  elastic 
symmetry.  The  C(a)  for  a  laminate  can  be  projected  from  its  mutually  perpendicular  planes  of 
elastic  symmetry  onto  the  composite  coordinate  system  (Xi,X2,X3)  about  the  X3  axis  via  the  fol¬ 
lowing  transformation  matrix. 


Q(a)  = 


C® 

0 

0 

0 

-2CS 


0 

0 

0 

ICS 


0 

0 

1 

0 

0 

0 


0 

0 

0 

c 

s 

0 


0 

0 

0 

-s 

c 

0 


cs 

-cs 

0 

0 

0 

(C^-S^) 


(13.8.2) 


where  C  =  cos  0,  S  =  sin  0. 

Then  the  transformed  constitutive  matrix  is  obtained  as 

C'(o)  =  Q(o)C(o)Q(a) 

In  general,  ^(a)  will  have  the  following  form 


(13.8.3) 
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(13.8.4) 


We  evaluate  the  C(Kp.eia„.),  i.e.,  the  constitutive  relation  for  the  macro  element  via  the 
expression 


^rep.elon.)  =  ™^(1)  +  (1~™)C'(2)  (13.8.5) 

where  m  =  %il%2  <  1- 

Now  C(jq,  eiem )  provides  us  an  effective  C  matrix  which  repeats  itself  in  the  composite  in 
its  thickness  direction.  In  other  words  it  represents  an  effective  material  in  which  the  effects  of 
various  fiber  orientations  and  different  material  properties  have  been  integrated  in  a  consistent 
manner. 

C(in*CTo)~  J  e]Q^)dX3  (13.8.6) 

0 


Fig.  13.8  Finite  element  mesh,  through  the  thickness 
Figure  13.8  shows  a  typical  representative  element,  composed  of  "n"  plies  (in  this  figure,  n  =  2). 

There  can  be  "m"  such  representative  elements  in  the  thickness  direction. 

Now  that  we  have  shown  a  consistent  way  of  arriving  at  the  macro  level  from  the  ply  level, 
we  are  in  a  position  where  we  can  talk  about  its  numerical  implementation.  There  are  a  number 
of  important  factors  which  have  to  be  considered  and  are  summarized  below.  These  technical 
issues  make  the  present  implementation  different  from  a  3D  continuum  model  and  also  from  the 
standard  2D  shell  elements,  whether  they  are  based  on  the  degenerated  shell  approach  or  the 
Cosserat  surface  approach. 

From  an  implementational  standpoint,  this  theory,  in  effect,  combines  the  2D  shell  effects 
with  the  3D  continuum  effects.  It  is  a  well  known  fact  that  a  shell  can  be  analyzed  as  a  congre¬ 
gate  of  3D  elements  in  which  the  in-plane  dimensions  of  the  finite  elements  are  of  the  order  of  its 
thickness  dimension.  The  drawback  lies  in  the  tremendous  storage  requirements  and  the  CPU 
intensive  calculations  of  the  large  systems  of  equations.  A  shell  formulation  allows  the  same 
engineering  accuracy  with  considerably  less  elements,  thereby  making  the  computations  cost 
effective  and  economic. 
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Following  is  a  brief  presentation  of  finite  element  formulation. 

13.9  Finite  Element  Formulation 

13.9.1  Strong  Form  of  the  Problem 

To  the  general  equilibrium  equations  we  add  the  boundary  conditions 


Uct  =  aa  on  (13.9.1a) 

U3  =  U3  on  (13.9.1b) 

6a  =  §o  o”  I'e  (13.9.1c) 

In  addition 

T^np=T‘onr“  (13.9.  Id) 

qX^lonPq  (13.9.  le) 

S“Pnp=§*  on  r,  (13.9.1f) 


where  corresponds  to  the  boundary  where  displacement  type  boundary  conditions  are 

applied;  F"  ^  Fq  corresponds  to  the  boundary  where  traction  type  conditions  are  applied;  Fg 
corresponds  to  the  boundary  with  prescribed  rotation;  F,  corresponds  to  the  portion  of  boundary 
with  applied  moments. 

As  usual 

n  =  noEa5no(Aa  (13.9.2) 

denotes  the  unit  outward  normal  to  the  boundary  dCl  of  domain  £2. 
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Multiplying  the  strong  form  of  the  problem  with  the  admissible  variations,  integrating  by 
parts  and  using  the  prescribed  essential  and  the  natural  boundary  conditions  we  arrive  at  the 
weak  form  of  the  problem. 

The  spaces  relevant  to  the  problem  are 


13-30 


where  denotes  the  space  of  square-integrable  functions  along  with  their  generalized 

derivatives  defined  over  fl,  and  H<|(£2)  is  the  subset  of  H*(Q)  whose  members  satisfy  zero  essen¬ 
tial  boundary  conditions. 

13.9.3  Finite  Element  Nomenclature 

Summary  of  Composite  Laminate  C*-Theory  Notation 


_  <Uj,U2,U3>^^ 

ew=ew=<0j,e2>('>r 

“  ®(Sp) 

T^2=ufi-8<^=ua-e« 

A3  0 


displacement  vector 
rotation  vector 

curvature  tensor 

shear  strain  vector 

moment  (stress  couple)  tensor 


x,“ 


0^^=  /  xSd^ 


u, 


a 


0« 


shear  force  vector 


prescribed  boundary  displacement 


prescribed  boundary  rotations 


The  variational  equation  for  the  Composite  Laminate  Theory  for  Shells  is  derived  from  the 
variational  equation  of  the  three-dimensional  theory  by  making  use  of  the  preceding  relations. 


Let  c  1^  be  a  bounded  open  set  with  piecewise  smooth  boundary  F;  n^  ^  1  denotes  the 
number  of  spatial  dimensions.  F  admits  the  following  decomposition 


FgUFh  =  F 


(13.9.5a) 
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r,nri,  =  <>  (13.9.5b) 

where  Fg  and  F],  are  the  pardons  of  the  boundary  with  prescribed  essential  and  natural  boundary 
conditions,  defined  as 

Fg  =  F„uFe  (13.9.6a) 

Fh  =  F,uFqUF,  (13.9.6b) 

1.  The  domain  £2  is  of  the  following  special  form 

n  =  {(x,y,z)  €  J?  I  z  €  [OjC|^] .  (x,y)  €  A  c  1^}  (13.9.7) 

2.  The  integrals  appearing  in  the  three  dimensional  variational  equations  are  replaced  by 

xr“ 

/••dfl=f  J  -  -dxjdA  (13.9.10) 

n  AO 

f  •••  dF  =  J<  *->dA+ J  j  -  dxjds  (13.9.11) 

U  A  dA 

where  ds  is  measured  along  the  perimeter  of  A. 

13.9.4  Finite  Element  Stiffness  Matrix  and  Load  Vector 

The  finite  element  stiffness  matrix  and  load  vector  may  be  obtained  directly  from  the 
matrix  form  of  the  variational  equation.  The  finite  element  approximations  for  u,u,0  and  0  are 
denoted  by  and  0**,  respectively. 

In  a  typical  element,  possessing  n^n  nodes, 

u‘‘=2Nt“i'  (13.9.12a) 

S=1 
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u‘=J  N.ii.‘ 

*=1 

(13.9.12b) 

®^=SN,e.*‘ 

t^\ 

(13.9.13a) 

^=IN.eJ 

(13.9.13b) 

«=1 


where  is  the  shape  function  associated  with  node  "a",  u. ,  ,  and  O]*  are  the  a*^  nodal 

values  of  u\  6*'  and  respectively.  It  is  not  necessary  to  assume  that  6^  and  be  defined 
in  terms  of  the  same  shape  functions  and  nodal  patterns.  However,  for  the  present  implementa¬ 
tion,  this  will  be  the  case. 

13.9,5  Nodal  Degrees  of  Freedom 


d*={dp*} 


d‘={d;} 


j 


“U 

n** 

U2, 

17** 

“3> 

ef. 


where  p  is  the  local  equation  number  at  the  element  level. 


(13.9.14a) 

(13.9.14b) 


(13.9.14c) 


(13.9.14d) 
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13.9.6  Matrix  Expressions 
Stress  Vectors 


The  resultant  stress  vectors  for  an  element  which  has  the  3D-cffccts  of  an  elasticity  element 
and  2D  effects  of  a  shell  like  element  are: 


x  = 


inplane 


4  = 


q.® 

qi® 


shear 


o  = 


off 

of? 

of? 


bending 


(13.9.15a) 


(13.9.15b) 


(13.9.15c) 


t=  {t^§)  through-the-thickness 


(13.9.15d) 


Strain  Vectors 

The  strain  vectors  corresponding  to  the  stress  vectors  are 


e  = 


rff  ■ 
10  • 
2Yf? 


inplane 


(13.9.16a) 


s= 


shear 


Kff 

K0  • 

2icf? 


bending 


(13.9.16b) 


(13.9.16c) 
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lyi?}  through-the-thickness  (13.9.16d) 

We  can  combine  the  strain  vector  for  inplane  effects  with  the  strain  vector  for  through-the- 
thickness  effects  to  yield  a  vector  that  incorporates  3D  effects 

e=- 

13.9.7  Matrix  Differential  Operators 

The  strain  vectors  can  be  written  in  terms  of  differential  operators  as  follows 


y(? 

Yg 

Y® 

Hi 


(13.9.17) 


aJ  a/ax* 
Aja/ax^ 
D'*’a/ax^ 

Afa/ax^+Aja/ax^ 


For  the  case  of  flat  geometry  and  linear  analysis 


(13.9.18) 


r  n 

01 

f  01 

oo 

II 

< 

—A _ 

[o\ 

,  D=- 

0 

1 1 J 

The  matrix  differential  form  for  the  modified  membrane  effect  then  becomes 


Bm  = 


a/ax^ 

a/ax^ 

a/ax^ 

a/ax^+a/ax^ 


Similariy,  for  the  shear  operator 


(13.9.19) 


So  =  7a3  =  Y  (“3  a  -  5o)  =  Y  (“3.0  -  6o) 
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Remark:  Note  the  introduction  of  the  rotation  field  in  the  above  strain-displacement  relation 


r  •N 


S„  =  [DTa/3X«  aJ] 


(13.9.20) 


Therefore  transverse  shear  strain  vector  S  is 


u 

e 


where 


D'^a/ax* 

D'^'a/ax^ 


for  flat  geometries: 


Bs  =  [Bsm  Bjb]  = 


a/ax^ 

a/ax^ 


Ar 

Aj 


(13.9.21) 


A 

Following  the  same  lines,  the  bending  strain  vector  S  which  comprises  of  the  operators 
(13.9.16c)  and  bending  membrane  coupling,  can  be  written  in  terms  of  matrix  diffoential  opera¬ 


tors: 


S-[Bbm 


f  ^ 


where 


Dja/ax^ 

A?"  a/ax^ 

Bbm  “ 

Dja/ax^ 

;  Bbb- 

A  J  a/ax^ 

Dja/ax^+Dja/ax^ 

A?  9/9x2 -i-Aja/ax^ 

(13.9.22) 


Specializing  for  flat  geometries 


13-36 


'ft 


=  0 


Bb„.  =  [0] 


a/3x^ 

dfdX^ 

d/ax^+d/ax* 


Finally,  a  total  matrix  differential  operator  B  can  be  defined  which  produces  the  total  strain  vec¬ 
tor  e  when  applied  to  the  displacement  field  u  and  the  rotation  field  6: 


e  =  B 


u 

6 


where 


B  = 


Bm  0 
Bsm  Bjb 

Bbm  Bbb 


(13.9.23) 


13.9.8  The  Strain  Displacement  Matrices 


We  can  write  the  matrix  differential  operators  in  tnms  of  element  shape  functions  as  fol¬ 
lows: 


Nm 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Nu 

Nu 

0 

0 

0 

0 

0 

-N. 

0 

0 

0 

Na 

0 

-N. 

0 

0 

0 

N..1 

0 

0 

0 

0 

0 

Ni^ 

0 

0 

0 

N..2 

Nu 

(1.3.9.24a) 


(13.9.24b) 


(13.9.24c) 


where  ”a"  stands  for  the  node  number. 
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Bm.  0 

K 

= 

Btm,  BjI^ 

K 

•  « 

Bbm,  Bbb^ 

13.9.9  Stiffness  Matrix 

Using  the  finite  element  assembly  operator,  tiie  stiffness  matrix  is  obtained  as 

K  =  ££K*  (13.9.25) 

hi  e::! 

where  L  is  tiie  total  number  of  finite  element  layers  through  the  thickness  and  is  the  total 
number  of  elements  in  each  layer. 

We  can  also  write  k  as 


K»A 


e=l 


nw 


=  A 


e=l 


(13.9.26) 


where  we  have  combined  the  inplane  membrane  effects  with  through  the  thickness  effects  to 
engender  the  modified  k^i,- 


Let  us  consider  each  of  the  stiffness  contributions  one  by  one. 
•  Modified  nembrane  stiffness 


5^^ 


where 
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'  L  " 

^^,=  J  /  BiC^B^jdCdn  (13.9.27) 

SxS  0  Sx4  4x4  4x5 

where  is  the  constitutive  matrix  for  membrane  effects  for  layer  /,  L  denotes  the  total  number 
of  finite  element  layers  and 


+  +1 

J  d  □  s  f  J  ’  *  *  d^dT)  Gamina  integral) 
a  -1-1 

j  is  the  determinant  of  the  Jacobian  defined  as 


j  =  dct 


*U  *i.C 
*2.11  *2,; 
*3.§  *3.n  *3.4 


We  can  write  the  above  equation  as 

=  (13.9.28) 

i=l 


is  the  integration  rule  required  to  exactly  evaluate  the  membrane  effects.  In  this 
equation  we  have  replaced  the  integral  sign  with  sumnaadon  over  all  the  integration  points  and 
have  also  multiplied  by  the  corresponding  weight  wH  Since  it  is  a  volume  integral,  so  the  loop 
over  the  nodes  covers  all  the  nodal  points  in  the  element  Furthermore  we  need  a  2  x  2  x  2 
integration  rule  to  evaluate  all  the  monomials  in  the  integral  exactly.  [Note:  2x2x2  rule  is  for 
an  8-node  brick  element]  We  need  to  use  the  appropriate  rule  for  higher  order  elements. 
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•  Bending  Stiffhess 


•^=n<22^i 

where  is  the  constitutive  matrix  for  the  bending  effects  for  layer  /,  and  j  is  the  determinant  of 
the  2D  reference  surface  Jacobian  defined  as 


j  =  dct 


*1.5  *1.11 

*2.5  *2.11 


t=l 

The  nomenclamre  used  in  the  jproceeding  expressions  is  defined  as  below. 

A  A 

fi,b;  finite  element  nodes  associated  with  the  reference  surface  (1  <=  S,  b  <  fien)- 

iiint((,,’  integration  rule  required  to  exactly  integrate  the  bending  contribution. 

B|,:  strain  displacement  matrix  associated  with  bending, 

fien-  number  of  nodes  associated  with  the  reference  surface. 


•  Shear  Stiffhess 


□ 

where  is  the  constitutive  matrix  for  the  shear  effects  for  layer  /.  Consequently 

C<'>  B^  J  w«  (13.9.30) 

isl 
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13.9.10  Evaluation  of  Stresses 


Let  us  consider  each  "layer  of  finite  elements"  through  the  thickness  direction.  Let  and 

represent  the  inverse  of  the  in-plane  and  out-of-plane  parts  of  the  compliance  matrices  for 
the  layers  /.  As  noted  previously,  we  are  proposing  a  "layerwise  shear  deformable  finite  element 
theory  for  thick  composites."  Consequendy,  the  bending  attribution  and  the  shear  contribution 
to  the  stiffness  matrix  are  obtained  by  summing  over  the  coireponding  attributions  for  the 
individual  finite  element  layers.  For  each  layer  /,  and  are  defined  as  follows: 


■  of?-.  ■ 

Kf? 

0^.  - 

.  =  ci')- 

Kg  ■ 

2«:S 

Cii  Ci2  Ci6 

Cl2  C22  C26 

C16  C26  Cgfi 


(0 


(13.9.31) 


Combining  the  in-plane  membrane  effect  with  the  through-the-thickness  effects  in  the  element 


f-inod(o  _ 
''memb  ~ 


of?- 


— 


«1.1 
«2.2 
U33 

«U+«2,1  J 


(13.9.32) 


where  is  the  modified  membrane  stress  for  finite  element  layer  /,  and  is  defined  as 


Cii  C12  Ci3  C16 
C12  C22  C23  C26 
Ci3  C23  C33  C35 
C16  C26  C36  Cgfi 


Similarly 
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f  y81 


(13.9.33) 


Q5 

C55 


(0 


13.9.11  External  Force  Vector 


Body  Force:  The  element  body  force  vector  is  given  as 


(body  ^ 

Sn«xl 

=  J  J  NjfjdCdO  (13.9.34) 

5x1  □  0 


where 


N.  0 
0  N, 
0  0 


0 

0 

N. 


-N.Z.ei,  N.Z.ei. 
-N.Z.e{„  N.Z.e;, 
~N.Z.e;^  N^,ei. 


Surface  Force:  The  element  surface  force  vector  is  defined  by 


fr^=/N7hj,X3da  ,  1X3  = 
□ 


X3"^  :  Top 
0  :  Bottom ' 


(13.9.35) 


where  X3  is  the  thickness  parameter,  h  is  the  surface  force  vector  (per  unit  surface  area)  and  j, 
defined  as 


J*  ~  ^ 


BASE 


13-42 


is  the  surface  Jacobian  of  the  reference  surface.  The  surface  force  includes  both  pressure  and 
shear  which  can  be  defined  as  below. 

a)  Pressure'.  In  our  case 

h  =  -Xjpn  (X3  =  0  or  Xj  =  X3”")  (1 3.9.36) 

lle^xe^ll 

where  p  is  the  pressure  and  n  is  the  unit  normal  vector  to  the  composite  surface. 

b)  Shear:  We  assume  that  the  shear  is  specified  in  the  ^  and  T[  directions  on  the  surface  in  ques¬ 
tion.  In  this  case  the  surface  force  vector  is  given  by 

hsh^e^-Hht^e,, 

where  h^  and  h^  are  the  shear  in  the  ^  and  ti  directions,  respectively. 

Edge  Force:  Suppose  we  wish  to  apply  a  distributed  loading  along  an  t] 
denote  the  distributed  surface  force.  The  nodal  forces  are 

+1X3“* 

f.*^**=J  J  (NThje)l(T,=H.ior-i)dCd^  (13.9.38) 

-1  0 

where 

je  =  1 1  X  I  I  (edge  surface  Jacobian) 

The  case  of  loading  along  an  ^  =  +1  or  -1  edge  is  handled  by  interchaing  ^  and  T|  in  the  above 
relations. 

Note  that  when  the  reference  surface  is  not  taken  to  be  the  midsurface,  nodal  moments  are 
produced  even  when  h  is  constant  (i.e.,  in  general,  f,^  *  0^^  *  0).  If  edge  forces  or  moments  are 
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=  +1  or  -1  edge.  Let  h 
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specified  per  unit  edge  lengdi,  then  nodal  forces  are  confuted  as  explained  next  Oxisider  an  T| 
=  +1  or  -1  edge.  Let  f/"  =  denote  the  edge  force  and  let  mi”  =  ini“(^)  denote  the  edge 
moment  The  nodal  forces  are  then  given  by 


-t-i 


-1 


fine 

M 

f|“ 

fine 


mi 

®2 


Ilx,^lld4 


(13.9.39) 


Note  diat  mj^  and  ml^  must  have  the  same  sense  as  0i  and  62.  The  result  is  made  applicable 
to  an  ^  =  +1  or  -1  edge  if  ^  and  ri  are  interchanged  in  the  above  relation. 


The  element  external  force  vector  is  thus  defined  by  adding  all  the  preceding  force  vectors 
as  follows 


pt  ^  jbody  jsurf  ^  (edge  (13.9.40) 

13.9.12  Boundary  Conditions 

It  is  important  to  realize  that  the  boundary  conditions  in  the  present  theory  are  not  always 
the  same  as  those  for  the  classical  thin  plate  theory.  The  differences  occur  in  the  specification  of 
the  "singly  supported"  case.  In  the  nxxlified  theory,  there  are  two  ways  of  going  about  this, 
depending  on  the  actual  physical  constraints.  See  Hughes  1987,  p.  324-327  fOT  the  necessary 
details. 
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14.0  NUMERICAL  SIMULATIONS 

14.1  Introduction 

The  finite  element  formulation  presented  in  the  previous  section  was  in^lemented  in  the 
ccMiiputer  program  FEAP  and  several  test  problems  were  analyzed  These  included  various  com¬ 
posite  laminatBS  with  flat  and  curved  geometries.  The  results  of  these  aiudyses  are  presented  in 
this  section.  For  each  simulation  the  geometry,  the  material  pr(^>erties,  the  boundary  c(»ditions, 
loadings,  and  the  finite  element  discretization  were  discussed  and  various  results  of  the  simula¬ 
tion  were  presented  at  the  end  of  this  section. 

Simulation  #1 

14. 1 . 1  Free  Edge  RnnnHarv.yalue  Problem  [45,-45]s 

The  first  numerical  simulation  is  a  prismatic  symmetric  laminate  having  traction-free  edges 
at  ;C-  i  a  and  surfaces  z  s  ±  h,  and  loaded  by  strain  applied  only  on  its  ends  at  y  »  constant. 
Each  layer  is  composed  of  unidirection  fiber-reinforced  material  such  that  the  fiber  direction  is 
defined  by  its  angle  0  with  the  y-axis. 

The  elastic  properties  of  various  composite  materials  used  in  this  numerical  simulation 
have  been  taken  from  N.  J.  Pagano  [p.  4,  Pagano  1989].  These  material  properties  are  E^^  = 
137.9  GPa,  Ej  =  E2  —  14.48  GPa,  Glx  =  Git ~ Grz  ~  5.86  GPa,  —  Vi^  =  ^tz  ~  0.21  where 
subscript  L  denotes  the  direction  parallel  to  the  fibers,  T  denotes  the  in-plane  direction  perpen¬ 
dicular  to  the  fibers,  and  the  subscript  Z  denotes  the  out-of-plane  direction. 

In  this  example  a  laminate  consisting  of  four  unidirectional  fibrous  composite  layers,  two 
with  their  axis  of  elastic  symmetry  (fiber  direction)  at  445  and  two  at  -45  to  the  longitudinal  lam¬ 
inate  axis  was  considered.  Figure  14.1  shows  the  lamiiuite  geometry  and  the  coordinate  system. 
1%  strain  in  opposite  directions  is  applied  aiy=0  and  y=L,  respectively,  while  it  is  restrained 
to  move  in  the  axial,  lateral  and  thickness  directions  at  the  plane  passing  through  y=  1V2.  In 
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order  to  solve  this  problem  a  finite  element  mesh  comprising  1920  composite  shell  elements  with 
2665  nodes  was  generated.  The  physical  dimensions  fOT  the  numerical  simulation  were  x  =  20,  y 
«  60,  z  s  2.5,  with  12  elements  in  ;cdirection,  40  elements  in  the  direction  and  4  finite  elements 
through  the  thickness.  For  each  finite  element  layer  through  the  thickness,  the  reference  surface 
was  assumed  to  be  associated  with  the  bottom  surface  of  that  layer.  In  the  following  figures  ele¬ 
ments  241-280  belong  to  the  bottom  layer  with  445°  ply  orientation,  elements  721-760  and 
1201-1240  belong  to  the  middle  two  layers  with  -45°  orientation  while  elements  1681-1720 
belong  to  the  top  layer  with  445°  orientation  at  a  section  cut  at  the  dashed  cross-section  in  Figure 
14.1.  In  these  results,  q-1  is  the  transverse  shear  stress  Oy,  transverse  shear  stress 

0X2-  The  stress  distribution  obtained  agrees  very  well  with  the  numerical  results  by  J.  N.  Pagano 
for  the  same  physical  dimensions  of  the  problem  and  material  properties  and  rather  much  refined 
spatial  discretization. 

1.  Figiues  1.1  and  1.2  show  the  various  stress  components  for  6  =  445°  and  0  s  -45°,  respec¬ 
tively.  For  a  deeper  insight,  please  see  Fig.  1.3-1.8  and  the  following  explanation. 

2.  Figure  1.4  shows  the  inplane  shear  stress  a,y  and  is  equal  in  magnitude  but  opposite  in 
sight  for  the  two  material  sets.  The  ratio  of  Oyy  and  o^y  agrees  closely  with  that  of 
Pagano’ s  numerical  simulation. 

3.  Figure  1.5  shows  the  through  the  thickness  stress  component  Gq.  It  shows  a  sudden 
increase  in  value  close  to  the  free  surface. 

4.  Figures  1.7  and  1.8  show  the  transverse  shear  stresses  and  Figures  1.9,  1.10  and  1.11 
represent  the  moments  generated  because  of  the  shearing  stresses,  thereby  giving  rise  to 
curvatures  near  the  firee  edges  of  the  laminate. 

5.  Figures  1.12,  1.13  and  1.14  represent  the  displacement  fields  in  the  axial,  through¬ 
thickness  and  the  transverse  directions  to  the  applied  loading,  respectively,  whereas  Fig. 
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Simulation  #2 

14.1.2  Free  Edee  Boundarv-value  Problem  with  a  circular  hole  [45,-45]s 

This  numerical  simulation  is  a  prismatic  symn^tric  laminate  with  a  circular  hole,  having 
traction-free  edges  at  ^  a  and  surfaces  z  ~  ±  h,  and  loaded  by  strain  ^)plied  only  on  its  ends 
at  constant  The  ratio  of  diameter  of  the  hole  to  the  width  of  the  composite  is  D/W  »  0.1. 
Each  layer  is  composed  of  unidirectional  fiber-reinforced  material  such  that  the  fiber  direction  is 
defined  by  its  angle  6  with  the  jf-axis.  The  matmal  properties,  and  the  physical  dimensions  of 
the  composite  are  the  same  as  in  the  preceding  case,  however  the  computational  mesh  is  quite 
different  To  solve  this  problem  a  mesh  containing  720  composite  elements  per  finite  element 
layer  and  with  four  layers  through  the  thickness  was  generated.  The  total  number  of  nodes  in 
this  problem  is  3780.  The  unit  diameter  cylindrical  hole  has  its  center  point  at  (0,0,0).  No  boun¬ 
dary  conditions  are  applied  on  the  surface  of  the  hole,  i.e.,  it  is  traction  free.  The  block  is  con¬ 
strained  to  move  in  the  axial,  transverse  or  latoal  direction  by  appropriately  constraining  the 
nodes  along  the  symmetry  lines  on  the  section  passing  through  y=0. 

1.  Figures  2.1  and  2.2  present  the  various  stress  con^nents  for  and  -45°,  respectively, 
at  a  7(z  plane  that  passes  through  ^  =  0. 

2.  Figure  2.3  shows  the  axial  stress  Oyy  across  the  width  of  the  laminate.  Elements  521-540 
belong  to  the  bottom  layer  with  445°  ply  orientation,  elements  1241-1260  and  1961-1980 
belong  to  the  middle  two  layers  with  -45°  orientation  while  elements  2681-2700  belong  to 
the  top  layCT  with  445°  orientation.  Results  are  shown  for  the  transverse  plane  along  the 
positive  ^direction  at  ^  =  0.  Away  from  the  hole  the  stress  distribution  obtained  agrees 
very  well  with  the  preceding  numerical  results  for  the  same  physical  dimensions  of  the 
problem  and  material  properties.  However  a  sharp  gradient  in  the  stresses  can  clearly  be 
seen  and  the  value  increases  to  almost  two  times  its  value  in  the  region  away  firom  the  hole 
or  the  free  edge. 
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3.  Figure  2.4  shows  the  inplane  shear  stress  a^y  and  is  eqiuil  in  magnitude  but  opposite  in  sign 
for  the  two  material  sets.  Once  again  close  to  the  hole  there  is  a  sudden  increase  in  the 
value  of  the  stress.  The  ratio  of  Cyy  and  o,y  agrees  closely  with  that  of  the  preceding 
numerical  simulation. 

4.  Figure  2.5  shows  the  through  the  thickness  stress  component  o^z.  Near  the  hole  the  value 
climbs  to  12  times  its  value  near  the  free  surface.  This  shows  that  if  the  composite  was 
designed  to  withhold  the  through  thickness  stress  of  the  free  edge  and  a  hole  was  drilled 
through  it.  the  very  high  stress  concentration  that  develops  around  the  hole  can  result  in 
delaminadon,  thereby  resulting  in  a  very  sharp 

5.  Figures  2.6  and  2.7  are  important  in  the  sense  that  they  show  the  various  stress  components 
along  the  circumference  of  the  circle  for  the  +45®  and  the  -45®  laminates,  respectively. 
These  stresses  have  been  evaluated  at  the  Gauss  points  that  are  closest  to  the  free  edge. 
The  details  for  the  major  stresses  can  be  seen  in  Hgs.  2.8-2.10. 

6.  Figures  2.11-2.14  represent  the  displacement  and  rotation  fields  along  the  radial  distance 
from  the  hole  in  the  lateral  direction  to  the  applied  strain.  They  should  be  compared  with 
Figs.  1.12-1.15  of  the  previous  numerical  simulation.  The  difference  between  the  two  sets 
of  plots  is  that  in  the  previous  case  we  had  plotted  the  entire  cross  section  while  in  the 
present  case  only  half  of  the  section,  i.c.,  radially  out  from  the  center  of  the  hole  is  shown. 
It  is  evident  that  in  the  region  away  from  the  hole,  the  numerical  results  agree  closely  with 
that  of  the  simulation  without  the  hole.  It  shows  that  the  edge  effects  that  arise  because  of 
the  circular  hole  are  localized  in  a  small  region  around  it  However,  the  important  observa¬ 
tion  is  that  although  this  region  is  small  relative  to  the  entire  domain,  the  gradients  in  these 
regions  are  the  steepest,  making  it  more  vulnerable  to  delamination  and  high  stress  concen¬ 
tration. 
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7. 


Figures  2.15-2.18  provide  a  deeper  insight  in  the  displacement  and  rotation  pattern  around 
the  circular  hole.  Figure  2.15  shows  the  axial  displacement  along  the  circumference  of  tht 
hole.  An  antisymmetric  deformation  can  clearly  be  seen  which  shows  that  an  analogous  3- 
D  defmmation  phenomenon  takes  place  at  the  free  edges  generated  by  the  hole. 


elem.  2681-2700 
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Simulation  #3 

14. 1.3  Rree  Edge  Boundarv-value  Problem  [0,90]s 

The  third  numerical  simulation  is  again  a  prismatic  symmetry  laminate  having  tiaction-fiee 
edges  at  ^  and  surfaces  2  =  ±h,  and  loaded  by  strain  applied  only  on  its  ends  at  y-  constant. 
Each  layer  is  composed  of  unidirectional  fiber-reinforced  material  such  that  the  fiber  direction  is 
defined  by  its  angle  6  with  the  j^axis. 

The  elastic  properties  of  various  composite  materials  used  in  this  numerical  simulation 
have  been  taken  from  N.  J.  Pagano.  In  this  example  a  laminate  consisting  of  four  unidirectional 
fibrous  composite  layers,  two  with  their  axes  of  elastic  symmetry  (fiber  direction)  at  0^  and  two 
at  90*^  to  the  longitudinal  laminate  axis  was  considered.  Figure  14.1  shows  the  laminate 
geometry  and  the  coordinate  system.  1%  strain  in  opposite  directions  is  applied  aty^O  and  y = 
L,  respectively,  while  it  is  restrained  to  move  in  the  axial,  lateral  and  thickness  directions  by 
appropriately  constraining  the  nodes  along  the  synunetry  lines  at  =  L/2  section.  The  physical 
dimensions  of  the  problem,  the  computational  mesh  and  the  boundary  conditions  are  the  same  as 
for  simulation  1,  and  therefore  are  being  omitted  here. 

1.  Figures  3.1  and  3.2  show  the  various  stress  components  for  6  =  0°  and  6  =  90°,  respec¬ 
tively.  For  a  deeper  insight  into  the  behavior  of  these  stresses  the  reader  is  referred  to  Figs. 
3.3-3.6. 

2.  Figures  3.7-3.10  represent  the  displacement  and  the  rotation  fields  developed  at  t/  =  L/2.  In 
Fig.  3.7,  die  antisymmetry  in  the  axial  displacement,  although  extremely  small,  can  still  be 
seen  which  shows  that  the  element  is  numerically  very  stable.  There  is  another  very  iiiqror- 
tant  observation  that  needs  to  be  made  at  this  stage.  Unlike  the  simulation  #1  with  [45,- 
4S}s  laminate  (Fig.  1.15)  in  which  the  edge  at  L  undergoes  a  twisting  deformation,  in 
the  [0,90]s  laminate,  undo*  the  same  boundary  and  loading  conditions  the  response  in  the 
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through  thickness  direction  still  stays  symmetric.  What  is  important  to  note  is  that  when¬ 
ever  one  is  working  with  ply  orientation  which  is  different  from  [0,90]s  case,  there  is 
always  going  to  be  a  3-D  phenomenon  at  the  edges.  So  for  all  practical  purposes,  in  com¬ 
posite  design  or  in  the  design  of  engineering  conqwnents  made  of  conqrosites,  if  the  load¬ 
ing  is  applied  at  a  certain  angle  to  the  ply  laminates  which  is  other  than  0  degrees  or  90 
degrees,  a  computational  tool  which  accounts  for  3-D  effects  is  a  must. 
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Simulation  #4 

14.1.4  Free  Edge  BnuTiHarv-value  Problem  with  a  Qroular  Hole  [0,90]s 

This  numerical  simulation  is  a  prismatic  symmetric  laminate  with  a  circular  hole,  having 
tracticm-free  edges  at  ;c-  ^  and  surfaces  z  s  ddi,  and  loaded  by  strain  applied  only  tm  its  ends  at 
y  -  constant  The  diameter  to  width  ratio  D/W  s  0.1.  The  material  properties,  and  the  physical 
dimensions  of  the  composite  are  the  same  as  in  the  preceding  three  cases.  The  computational 
mesh  is  identical  to  the  one  used  in  Simulation  2.  The  difference  is  that  in  the  present  case  it  is 
[0,90]s  laminate.  The  unit  diameter  cylindrical  hole  has  its  center  point  at  (0,0,0),  see  Fig.  14.2. 
No  boundary  conditions  are  applied  on  the  surface  of  the  hole,  i.e.,  it  is  traction  free.  The  block 
i  >  constrained  to  move  in  the  axial,  transverse  at  lateral  direction  by  appropriately  constraining 
the  nodes  along  the  syimnetry  lines  at  y^O  surface.  We  have  modelled  the  entire  composite 
without  the  assumption  of  symmetry. 

1.  Figures  4.1  and  4.2  show  the  relative  magnitudes  of  all  the  stress  components  for  the  layers 
at  0  degrees  and  90  degrees,  respectively.  As  can  be  seen,  the  axial  stress  is  orders  of  mag¬ 
nitude  greater  than  all  the  other  stresses  and  close  to  the  hole  it  increases  to  alnwst  three 
times  its  value  than  in  the  rest  of  the  domain. 

2.  Figures  4.3  and  4.4  represent  the  axial  stress  Oyy  and  inteimalinar  normal  stress  at  the 
periphery  of  the  circular  hole. 

4.  Figures  4.5-4.6  and  Figs.  4.7-4.9  represent  the  displacement  field  along  the  radial  distance 
from  the  hole  in  the  transverse  direction  to  the  applied  loading,  and  along  the  circumfer¬ 
ence  of  the  circular  hole,  respectively. 
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14.1.5  Extension  Test:  Comparison  with  J.  N.  Reddv  [45,*45]s 

We  have  selected  this  numerical  from  the  paper  of  J.  N.  Reddy.  The  same  problem  with 
identical  boundary  conditions,  material  properties  and  ply  (mentations  has  also  been  solved  by 
R.  M.  Jones  in  his  book  on  Mechanics  of  Composite  Materials. 

Consider  a  thick,  symmetric,  angle-ply  laminate  [45,-45]s  subjected  to  axial  displacements 
on  the  ends.  The  laminate  has  a  length  of  2L,  width  2W,  and  thickness  2h,  with  L  =  lOW  and  W 
=  8h  (see  Fig.  14.1).  Each  of  the  four  material  layers  is  of  equal  thickness  h,  and  is  idealized  as  a 
homogeneous  orthotropic  material  with  the  following  properties  expressed  in  the  material  coor¬ 
dinate  system: 


EL*20xl0^psi  .  Br  =  Ej=:2.1xl0^psi 
Glx  ~  Glz  “  ®ZT  ~  0.85  X  10^  psi 
l^LT  ~  l^LZ  IZ-  0«21 

where  subscript  L  denotes  the  direction  parallel  m  the  fibers,  subscript  T  denotes  the  inplane 
direction  perpendicular  to  the  fibers,  and  the  subscript  z  denotes  the  out-of-plane  direction.  The 
origin  of  the  global  ccxndinate  system  coincides  with  the  centroid  of  the  three-dimensional  com¬ 
posite  laminate.  The  ;^coordinate  is  taken  along  the  width  of  the  laminate;  the  ^(xxndinate  is 
taken  along  the  length  of  the  laminate;  and  the  z-coordinate  is  taken  through  the  thickness  of  the 
laminate.  Since  the  laminate  is  symmetric  about  die  3(y-planc,  only  the  upper  half  of  the  lam¬ 
inate  is  nxxlelled  Thus  the  computational  domain  is  defined  by  (-W  ^  W,  -L  ^  ^  L,  0  ^  z 
^  2h).  The  displacement  boundary  conditions  for  this  problem  are 


Ui(0,-L,0)  =  0  Ui(0JL,0)  =  0 
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U2(;&-L,z)  *  0  U2(^z)  =  Uo 

U3(%y0)  =  0 

In  order  to  solve  this  problem  a  finite  element  mesh  comprising  1920  coixq)Osite  shell  ele¬ 
ments  with  2665  nodes  was  generated.  The  physical  dimensions  for  the  numerical  simulation 
were  L  s  200,  W  «  20,  z  s  2.5,  with  40  elements  in  redirection,  12  elements  in  the  direction 
and  4  ctmiposite  elements  through  the  diickness. 

1.  Figures  5.1-5.3  represent  the  various  significant  stress  components  through  the  cmnposite 
cross  section  cut  at  ^  =  0,  i.e.,  at  midspan  of  the  loading  axis.  These  numerical  results  are 
very  close  to  the  computed  numerical  values  of  I.  N.  Reddy  and  R.  M.  Jones,  although  both 
are  using  a  much  higher  refinement  of  the  mesh  in  the  thickness  direction  when  compared 
to  our  spatial  discretization. 

2.  Figures  5.4-5.6  represent  the  bending  stresses.  The  general  behavior  is  the  same  as  of  the 
previous  cases  for  [45,-45]s. 

3.  Figures  5.7-5.9  represent  the  displacement  fields  as  obtained  at  y  =  0  and  y  =  L,  respec¬ 
tively.  It  is  interesting  to  note  that  the  twisting  in  the  through-thickness  displacement  field 
that  we  observed  in  Simulation  1  is  also  present  in  the  present  simulation. 
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Simulation  #6 

14.1.6  Bending  Analysis  [45,-45, 45]s 

This  simulation  presents  a  bending  test  involving  a  square  plate  with  uniformly  distributed 
load.  In  diis  simulation  the  complete  plate  without  the  assumption  of  symmetry  was  analyzed. 
The  plate  is  composed  of  three  laminates  stacked  in  45/-45/45  cross-ply  construction  with  6 
measured  from  the  ^  axis.  The  material  properties  of  Simulation  1  which  were  taken  firom  J.  N. 
Pagano  have  been  used  in  this  case  as  well.  The  finite  element  mesh  is  composed  of  675  ele¬ 
ments  with  225  elements  per  finite  element  layer  and  3  finite  element  layers  through  the  thick¬ 
ness.  The  distribution  of  elements  is  15  x  15  x  3.  The  total  number  of  nodes  is  1024. 

In  the  present  numerical  simulation  the  plate  is  clamped  along  the  boundaries  and  a  uni¬ 
formly  distributed  load  of  intensity  q,,  is  applied  in  the  +z  direction.  Numerical  results  are 
reported  for  the  section  cut  &ttf=  b/2. 

1.  Figures  6.1  and  6.2  present  the  various  components  of  stresses  for  the  445  ply  on  the  ten¬ 
sion  side  and  the  -45  ply  which  lies  on  the  axis  of  symmetry  in  the  z  direction,  respectively. 
A  detailed  analysis  of  stresses  can  be  seen  in  Hgs.  6.3-6.8.  Since  the  edges  are  clamped,  the 
stresses  have  a  certain  finite  value  at  the  boundary. 

2.  Figures  6.9  and  6.10  present  the  bending  moments. 

3.  We  present  the  displacement  and  the  rotation  fields  in  Figs.  6.1 1-6.13. 

We  would  like  to  mention  that  bending  analysis  similar  to  this  simulation  has  not  been  reported 
in  the  literature. 

Figure  6.3  shows  bending  stresses  in  x  direction  (o„)  for  top,  center  and  bottom  layers  of 
the  plate.  For  the  center  layer,  the  stresses  are  calculated  at  the  Gaussian  integration  points 
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which  are  located  at  the  neutral  axis  of  the  plate  and  the  finite  element  simulations  predict  the 
correct  zero  values.  Stresses  for  the  top  and  bottom  layers  are  symmetric  with  regard  to  the  neu¬ 
tral  axis  but  witii  opposite  sign. 

Bgure  6.4  presents  bending  stresses  in  y-y  directions  for  the  same  cross-section. 

Figure  6.S  is  the  distribution  of  in-plane  shear  for  the  three  plies. 

Figure  6.6  presents  the  distribution  of  transverse  normal  stress  component. 

Figures  6.7  and  6.8  show  the  shear  stress  components  for  various  layers. 
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Figure  14.3 
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Sinuilation  #7 

14.1.7  Bending  Analysis  [45,-45,45]s 

This  simulation  again  presents  a  bending  tost  involving  a  square  plate  with  uniformly  distri¬ 
buted  load,  but  with  Simply  Supported  edges.  Once  again  in  this  simulation  the  complete  plate 
without  the  assumption  of  symmetry  was  analyzed.  The  plate  is  composed  of  three  laminates 
stacked  in  45/-45/45  cross-ply  construction  with  0  measured  fimn  the  axis.  The  physical 
dimensions,  material  properties,  loading  conditions  and  the  finite  element  mesh  used  for  analysis 
remain  unchanged  firom  simulation  #6. 

1.  A  detailed  analysis  of  stresses  can  be  seen  in  Figs.  7. 1-7.6.  Unlike  the  clamped  case,  in  the 
singly  supported  case  the  stresses  tend  to  vanish  at  the  free  edges. 

2.  Figures  7.7-7.9  present  the  bending  stresses. 

3.  We  present  the  displacement  and  the  rotation  fields  in  the  remaining  plots. 
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Simulation  #8 

14. 1.8  Cylindrical  Shell  with  Free  Edae  Boundary  [4S-4S]s 

The  last  numerical  simulation  is  a  cylindrical  shell  having  traction-£ree  edges  at  r  «  ri^2> 
loaded  by  1%  strain  applied  in  -z  direction  on  its  ends  at  9  s  90  degrees.  Each  layer  is  composed 
of  unidirectional  fiber-reinforced  material  such  that  the  fiber  direction  is  defined  by  its  angle  0 
with  the  y-axis. 

The  elastic  properties  of  various  composit  materials  used  in  this  numerical  simulation  have 
been  taken  from  N.  J.  Pagano  [p.  4.  Pagano  1989].  In  this  example  a  laminate  consisting  of  diree 
unidirectional  fiberous  composite  layers,  the  top  and  bottom  with  their  axis  of  elastic  symmetry 
(fiber  direction)  at  445  while  the  center  one  at  -45  to  the  longitudinal  laminate  was  considered. 
Figure  14.4  shows  the  laminate  geometry  and  the  coordinate  system.  1%  strain  in  the  -z  direc¬ 
tions  is  applied  at  the  cross-section  at  9  90,  while  it  is  restrained  to  move  along  the  cross  sec¬ 
tion  at  9  s  0.  In  order  to  solve  this  problem  a  finite  element  mesh  comprising  675  composite 
shell  elements  with  1024  nodes  was  generated.  The  physical  dimensions  for  the  numerical  simu¬ 
lation  were  r^  =  20,  r2  =  21.5,  z  =  20,  with  15  elements  in  the  9  direction,  15  elements  in  the  y 
direction  and  3  composite  elements  through  the  thickness.  For  each  finite  element  layer  through 
the  thickness,  the  reference  surface  was  assumed  to  be  associated  with  the  bottom  surface  of  the 
composite  shell.  We  will  be  looking  at  a  section  cut  at  Y  =  0  for  the  displacement  and  the  stress 
fields. 

1.  Figure  8.1  shows  across  the  width  of  the  laminate.  Elements  1-15  belong  to  the  bottom 
layer  with  445  ply  orientation,  elements  226-240  belong  to  the  middle  layer  with  -45  orien¬ 
tation  while  elements  451-465  belong  to  the  top  layer  with  445  orientation. 

2.  Figure  8.2  shows  Oyy.  Figure  8.3  shows  stress  component  o^.  At  9  =  0  it  shows  tension  on 
the  inner  layer  and  compression  on  the  outer  layer  with  445  ply  orientation. 
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3.  Figures  8.4-8.6  represent  the  above  three  stress  components  at  a  s^cion  cut  at  y  =  b/2. 

4.  Figures  8.6-8.8  represent  the  displacement  field  at  section  y  =  0  and  Figs.  8.9-8. 1 1  represent 
the  corresponding  values  at  an  interior  section  with  y  =  b/2. 


Figure  14.4 
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Stresses, [45, -45]s,1%  strain  along  Y 
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Figure  8.2 
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Figure  8.3 


Sig-11,[45,-45,45]  Cylindrical  Shell 
Section  cut  at  y  =  b/2 

40000 
30000 
20000 
10000 
0 

-10000 
-20000 
-30000 
-40000 

0  10  20  30  40  50  60  70  80  90 

Distance  along  Circumference  (Deg.  Phi) 


el  106-120, [45]  ♦  el  331-345,[-45]  el  556-570,145] 

Figure  8.4 


Stress  (KPa)  Stress  (KP^ 


14-66 


Sig-22.[45.-45.45]  Cylindrical  Shell 
Section  cutat  Y  =  b/2 
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Figure  8.5 


Sig-33, [45,-45,45]  Cylindrical  Shell 
Section  cut  at  y  =  b/2 
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Figure  8.6 
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Figure  8.7 


Uy.[45, -45.45],  Cylindrical  Shell 
Section  cut  at  Y  =  0  (free  edge) 
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Figure  8.8 
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Uz.  [45.-45.451,  Cylindrical  Shell 
Section  cut  at  Y  =  0  (free  edge) 
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Figure  8.9 


Ux, [45,-45,45],  Cylindrical  Shell 
Section  cut  at  Y  =  b/2 
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Figure  8.10 
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